Cite this Article
Hu Gan, Huang Jian-Qi, Wang Ya-Ning, Yang Teng, Dong Bao-Juan, Wang Ji-Zhang, Zhao Bo, Ali Sajjad, Zhang Zhi-Dong. Unconventional lattice dynamics in few-layer h-BN and indium iodide crystals. Chinese Physics B, 2018, 27(8): 086301  
Permissions
Unconventional lattice dynamics in few-layer h-BN and indium iodide crystals
Hu Gan1, 2, Huang Jian-Qi1, 2, Wang Ya-Ning1, 2, Yang Teng1, 2, †, Dong Bao-Juan1, 2, Wang Ji-Zhang1, 2, Zhao Bo3, Ali Sajjad1, 2, Zhang Zhi-Dong1, 2
Shenyang National Laboratory for Materials Science, Institute of Metal Research, Chinese Academy of Sciences, Shenyang 110016, China
School of Materials Science and Engineering, University of Science and Technology of China, Hefei 230026, China
College of Sciences, Liaoning Shihua University, Fushun 113001, China

 

† Corresponding author. E-mail: yangteng@imr.ac.cn yanghaiteng@msn.com

Project supported by the Major Program of Aerospace Advanced Manufacturing Technology Research Foundation from NSFC and CASC, China (Grant No. U1537204), the National Key Research and Development Program of China (Grant No. 2017YFA0206301), and the National Natural Science Foundation of China (Grant No. 51702146).

Abstract

Unusual quadratic dispersion of flexural vibrational mode and red-shift of Raman shift of in-plane mode with increasing layer-number are quite common and interesting in low-dimensional materials, but their physical origins still remain open questions. Combining ab initio density functional theory calculations with the empirical force-constant model, we study the lattice dynamics of two typical two-dimensional (2D) systems, few-layer h-BN and indium iodide (InI). We found that the unusual quadratic dispersion of flexural mode frequency on wave vector may be comprehended based on the competition between atomic interactions of different neighbors. Long-range interaction plays an essential role in determining the dynamic stability of the 2D systems. The frequency red-shift of in-plane Raman-active mode from monolayer to bulk arises mainly from the reduced long-range interaction due to the increasing screening effect.

PACS: 63.22.-m;63.20.D-;63.22.Np
Keyword:flexural mode;quadratic dispersion;long-range interaction;frequency redshift;dielectric screening
1. Introduction

Two-dimensional crystals have been attracting tremendous attention in the last few decades due to its great potential in fundamental research and multitude of applications. The quantum confinement due to dimension reduction gives rise to rich and novel physical properties which cannot be obtainable in their three-dimensional (3D) counterparts. It is very essential to pinpoint and understand the novel properties before it can be fully explored for diverse applications. So far there have been some phenomena which have not been fully understood and still await for a definitive answer, for example, the unusual quadratic phonon dispersion of flexural mode near the Brillouine zone center,[1–7] the counterintuitive frequency red-shift of in-plane mode with increasing layer thickness.[8–11]

The quadratic phonon dispersion behavior (ω ∼ q2) of the flexural (out-of-plane) acoustic mode against the linear dispersion (ω ∼ q) of the other acoustic phonon modes in 2D crystals have been observed widely ever since debut of graphene.[1] Due to the lack of a good understanding on the novel dispersion relation, there is a controversy on the contribution of the flexural mode to lattice thermal conductivity.[2–5] Liu et al.[6] studied such a quadratic dispersion relation of quasi-two-dimensional structures based on continuum elasticity theory, but the microscopic mechanism is still unknown. Jiang et al. studied microscopically the effect of rotation symmetry on the nonlinear dispersion of the flexural mode.[1,7] However, this analysis may be only limited in graphene system, not applicable to other 2D systems.

Another issue we are concerned with here is the unconventional frequency red-shift of in-plane Raman mode with increasing number of layers. It was firstly observed experimentally in MoS2 by Lee et al.[8] and also in h-BN crystals by Gorbachev et al.[9] Molina–Sanchez et al.[10] theoretically calculated the long-range screening as a function of layer thickness of MoS2 and WS2 crystals, and ascribed the thickness-reduced frequency of the in-plane mode to the enhanced screening effect with increasing layer number. However, Song et al.[11] considered that so-called Davydov splitting will lift symmetry and some degenerate mode will be split into multiple modes with different frequencies. Obviously the physical origin for this frequency red-shift is still unclear.

This study will focus on the lattice dynamics of two typical 2D systems, few-layer h-BN and indium iodide (InI), to get an insight of the two aforementioned phenomena by combining ab initio density functional calculations with the force-constant model. We found that the unusual quadratic dispersion of vibrational frequency on wave vector may be comprehended based on the competition between atomic interactions of different neighbors. Long-range interaction plays an essential role in determining the dynamic stability of the 2D systems. The frequency red-shift of in-plane Raman-active mode from monolayer to bulk is mainly due to the reduced long-range interaction by the increasing screening effect.

2. Methods

We used ab initio density functional theory as implemented in the Quantum Espresso package.[12] A periodic boundary condition with monolayer structures represented by a periodic array of slabs separated by a vacuum region(∼ 15 Å). We used the Perdew–Zunger exchange–correlation functional[13] within the local-density approximation (LDA) and norm-conserving pseudopotentials.[14] The Brillouin zone of the primitive unit cell of the 2D structures was sampled by 13 × 13 × 1 k-points.[15] The kinetic energy cutoff for wave functions was set at 150 Ry (1 Ry = 13.6056923 eV) and 10−10 Ry for a total energy difference between subsequent self-consistency iterations as the criterion for reaching self-consistency. All geometries are optimized using the conjugate gradient method,[16] until none of the residual Hellmann–Feynman forces exceeds 4 × 10−5 Ry/Bohr.

The phonon dispersion relation of InI and h-BN was calculated based on two methods, one is density functional perturbation theory (DFPT)[17] and the other the force constant model. The non-resonant Raman intensities were calculated based on Placzek approximation as introduced by Lazzeri and Mauri.[18]

The inter-atomic long-range interaction of polar materials includes the non-analytic term in the lattice dynamical matrices, which comes from the dipole–dipole interactions as an effect of vibration-induced polarization. Such a non-analytical behavior is available through the knowledge of the Born effective charges tensor Z* and the electronic dielectric permittivity tensor ɛ∞. The non-analytical, direction-dependent term at the vicinity of BZ center can be expressed by the following equation:[19,20]

The following ansatz of the dipole–dipole force constants or long-range force constants can well reproduce the aforementioned non-analytical behavior:[19,20]
where Δα = Σβ (ɛ−1)αβdβ is the conjugate of the vector d relating nuclei, while the norm of the latter is . The contribution of these long-range force constants to the dynamical matrix can be calculated using the Ewald summation technique.[19,20]

The force constant model is also used to calculate the phonon dispersion relation. Within this model, interatomic interactions up to the 4th neighbor are considered. The force constant tensor is defined as follows:

in which , , and are the in-plane longitudinal and transverse force constant and the out-of-plane transverse force constant, between the n-th atom in the unit cell and the n′-th atom at its m-th neighbor, under the in-plane longitudinal (ɛαɛβ), transverse displacement ( ), and the out-of-plane transverse displacement ( ), respectively.

3. Results and discussions

The atomic structures of monolayer and bulk indium iodide (InI) and hexagonal BN (h-BN) are shown in top view and side view, respectively, in Fig. 1. Different from atomic thin monolayer h-BN, monolayer InI has two sub-layers with In and I interchanged. The structures of h-BN and InI few layers are also different. For h-BN, AA’ stacking is energetically preferred with B(N) of one layer on top of N(B) of the other, while in bulk InI, Indium atom in one layer sits on top of the square center of the other layer due to repulsion between lone pairs of indium.[21–23] Two systems with different monolayer structures and stacking order for few layers may be representative to study long-range interaction effect on their lattice dynamics.

Fig. 1. (color online) Atomic structure of (a) InI and (b) h-BN in the monolayer and bulk form.
3.1. Lattice dynamics of h-BN

To study their lattice dynamics, it is essential to calculate the phonon dispersion relation. Figure 2 shows the phonon dispersion relation of both monolayer and bulk h-BN and the layer number N-dependence of Raman shifts. As seen from Fig. 2(a), there exists an out-of-plane acoustic phonon mode (ZA) (highlighted in red) with quadratic dispersion in both monolayer and bulk forms of h-BN. Raman-active mode (∼ 1380 cm−1) in zone center, as marked with blue filled dot in Fig. 2(a), is a bond-stretching mode, as shown in Fig. 2(b). Along with the increasing number of layers, the vibration mode is split into multiple modes due to degeneracy lifting by interlayer interaction (Davydov splitting). Interestingly, all frequencies go redshift.

Fig. 2. (color online) (a) Phonon dispersion relation of monolayer and bulk h-BN crystal. Black lines are from DFPT calculation, and red line of the ZA mode is a fit by the force constant model. (b) E2g mode and its splitting as a function of layer number N. The frequency value is written in bracket.

To have a better view of such an unusual behavior, we show in Fig. 3 the N-dependence of Raman shift and intensity of the E2g mode of h-BN. Not only Raman frequency, but also Raman intensity, changes with number of stacking layers of h-BN. Dashed red line shows the thickness-induced redshift of Raman shift, as well as the increasing Raman intensity. In the inset of Fig. 3, we compare the results from the experiment[9] (blue dots) and our calculation (red dots). Note that the intensity is normalized by the intensity value of 4-layer structure, and Δω is defined as the frequency difference between N-layer ω(N) (N = 1, 2, 3, . . . ) and bulk ω(∞). A good agreement is reached between experiment and calculation.

Fig. 3. (color online) Calculated Raman shift and intensity of the E2g mode as a function of layer number N of h-BN crystal. The inset shows the layer number dependence of both the peak intensity and frequency change Δω. Δω(N) is defined as the frequency difference between N-layer ω(N) and bulk ω(∞). Raman intensity I(N) is normalized by the Raman intensity at N = 4. Experimental data in blue filled circle are compared with the calculated ones in red filled square. The line-width of 12.00 cm−1 is used.
Table 1.

Screening effect of h-BN as a function of layer number. ψlr and ψsr are long-range and short-range force constants, respectively. Z* is Born effective charge.

.

The linear N dependence of Raman intensity observed in many 2D materials can be well understood from the quantized optical absorption[24,25] (∼ Nπα with α = e2/hc ≈ 1/137, N as the number of layers). The thickness-induced frequency redshift of in-plane mode has also been observed in other 2D systems such as MoS2, WS2, and MoTe2.[8,10,11] It is still in argument and not clear whether this is due to enhanced screening effect on the long-range interaction[8,10] or by Davydov splitting.[11] From our calculation results, we found that, unlike the split frequency values which go both up and down with increasing layer number in MoTe2, the split frequency values all go down, which excludes the origin from the Davydov splitting.

To unveil the possible origin due to the screening effect, we obtained both long-range and short-range lattice force constants, as well as dielectric constant as a function of layer number N of h-BN from ab initio calculations, as shown in Table 1. From Table 1, with increasing layer number N, the dielectric constant goes up, suggesting an increased screening effect on both in-plane and out-of-plane directions due to more layers. Meanwhile, the total force constant, which consists of both long-range and short-range interactions, decreases. In contrast to TMDCs[10] in which short-range interaction increases with N while long-range interaction decreases with N, both the long-range and short-range parts in h-BN change in the same way with N, indicating comparable role played by both interactions.

Now we have a look at the possible origin behind the quadratic dispersion of out-of-plane acoustic (ZA) mode in h-BN. In Fig. 4, we show the dispersion relation of ZA mode of h-BN as a function of out-of-plane force constant due to the 2nd neighbors. As discussed in details in Appendix A, to guarantee that the quadratic dispersion of ZA, an opposite sign has to remain between of the nearest neighbor and of the 2nd nearest neighbor, and . Here in h-BN case, dyn/cm and dyn/cm (1 dyn = 10−5 N), which meets the condition. Clearly, the dispersion changes from a linear to quadratic behavior when varying from positive (same sign as the ) to negative (opposite sign from the ). Basically the competing relation between the different orders of neighbors is essential for an appearance of the quadratic dispersion of ZA mode.

Fig. 4. (color online) (a) Definition of all the components of force constant. with the subscript a = r, ti, to represent in-plane radial, in-plane and out-of-plane transverse, respectively, n represent the n-th nearest neighbor (nn) and the superscript A, B represent the element B at the n-th nn of the element A. Color circles are used to show the n-th nn. (b) Dispersion relation of out-of-plane ZA mode of h-BN monolayer as a function of out-of-plane force constant or due to the 2nd nearest neighbors. Dispersion from a linear to quadratic behavior is obtained as varying or while keeping other parameters fixed (i.e., ). The unit of force constant is 104 dyn/cm (1 dyn = 10−5 N).
3.2. Lattice dynamics of InI

Monolayer InI was theoretically predicted to be stable by Wang et al.[23] It has two atomic layers of alternating indium and iodine atoms with the In atom extruding slightly (≈ 0.50 Å) out of the layer plane. Few-layer and bulk phase have orthorhombic structure of Cmcm symmetry (D2h point group). Here we studied bi-layer, tri-layer, and bulk besides monolayer structure. Phonon dispersion relation for both monolayer and bulk is shown in Fig. 5(a). The anisotropy of phonon dispersion relation along X and Y directions in bulk InI is due to the anisotropic interlayer stacking configuration, in which Indium atoms form a zigzag chain between two sub-layers along one direction, as seen more clearly in Fig. 1(a). Out-of-plane acoustic ZA mode is highlighted in red. A quadratic dependence on wave vector is observed in monolayer while a linear behavior in bulk InI. Raman active modes are denoted by the irreducible representation Bg and Ag, which represent in-plane and out-of-plane modes, respectively.

Fig. 5. (color online) (a) Phonon dispersion relation of monolayer and bulk InI crystal. BZ with the high-symmetry points is shown. (b) Visualization of typical in-plane Bg mode and out-of-plane Ag mode of monolayer, bilayer, trilayer and bulk InI.

Two typical Raman active vibrational modes are studied here, one is the layer breathing Ag mode having the highest frequency value and the other the shearing Bg mode having the lowest frequency value, as highlighted in green and blue dots, respectively, in Fig. 5(a) and also visualized in Fig. 5(b). Since monolayer consists of two sublayers, it is reasonable to have both layer breathing and shearing modes starting from monolayer to few layers.

The N dependence of both Ag and Bg modes of InI is studied in Fig. 6. For the out-of-plane layer breathing Ag mode, it has relatively large Raman intensity, as seen from Fig. 6(a). Both Raman shift ω and intensity increase with increasing layer number N, which is summarized in Fig. 6(b). Similar phenomena have been seen in layered structures with unit cell thicker than one atomic layer, such as TMDCs.[8,10] Interestingly, the curve of ω versus inverse N can be fitted well by a function of (with A = 92.15 cm−1 and B = 0.578, as plotted by a black dashed line in Fig. 6(b)), which is reminiscent of the linear-chain model [26,27] and suggests a short-range interaction. Note that “B = π” appears in the inter-layer mode ( ) in h-BN and graphene, while in InI, B ≠ π for the layer breathing mode. We may ascribe this behavior to the unit-cell structure of InI which has two sub-layers other than one single atomic layer in h-BN and graphene. This structural difference may invalidate the simple chain model.

Fig. 6. (color online) (a) and (c) Calculated N-dependent Raman spectra. The line-width used is 0.20 cm−1. Green dashed frame in panel (a) and dashed arrow in panel (c) highlight the typical out-of-plane Ag mode and in-plane Bg mode, respectively. Raman shift and intensity of (b) Ag and (d) Bg modes as a function of inverse layer number 1/N of InI crystal. A black dashed curve in panel (b) is a fit to ω(1/N) data. To note that logarithmic scale is used for intensity in panels (b) an (d).
Table 2.

Screening effect of InI as a function of layer number. ψlr and ψsr are long-range and short-range force constants, respectively. Z* is Born effective charge.

.

In contrast, the N dependence of the lowest Bg shearing mode in Figs. 6(c) and 6(d) shows that Raman shift goes down while Raman intensity goes up with the increasing N, similar to that of E2g mode in h-BN crystals and TMDCs,[8,10] suggesting a suppressed long-range interaction due to thickness-enhanced dielectric screening effect. This can be seen more clearly from Table 2.

We further studied the source of quadratic dispersion of ZA mode frequency of InI monolayer. As shown in Fig. 7, a conversion from linear to a quadratic dispersion is realized once a sufficiently large force constant due to the 2nd neighbors competes with next nearest interaction with a negative sign. In the meantime, a concave shape of dispersion at the edge corner of Brillouine zone appears, which is coherent with what is obtained from ab initio DFT calculation, as seen from Fig. 5(a).

Fig. 7. (color online) (a) Definition of all the components of force constant. (b) Dispersion relation of ZA mode of InI as a function of force constant or at the 2nd nn. Dispersion from a linear to quadratic behavior is obtained with varying or .
4. Conclusion

In summary, we have studied the N dependent lattice dynamical properties of indium iodide (InI) and h-BN crystals, by combining ab initio density functional theory calculations with the empirical force-constant model. Unconventional thickness-induced frequency red-shift of in-plane vibrational modes (E2g in h-BN and Bg in InI) has been found from our calculations, and the theoretical results agree quite well with the experimental results available. This behavior can be well comprehended based on the dielectric screening effect, which is enhanced with increasing layer number and therefore can suppress long-range interaction responsible for those in-plane Raman active modes. Meanwhile, unconventional quadratic dispersion of the out-of-plane flexural ZA modes appears due to the sufficient competition between interaction among different nearest neighbors.

Reference
[1] Jishi R Venkataraman L Dresselhaus M Dresselhaus G 1993 Chem. Phys. Lett. 209 77
[2] Mounet N Marzari N 2005 Phys. Rev. 71 205214
[3] Lindsay L Broido D A Mingo N 2010 Phys. Rev. 82 115427
[4] Gunst T Kaasbjerg K Brandbyge M 2017 Phys. Rev. Lett. 118 046601
[5] Nika D L Ghosh S Pokatilov E P Balandin A A 2009 Appl. Phys. Lett. 94 203103
[6] Liu D Every A G Tománek D 2016 Phys. Rev. 94 165432
[7] Jiang J W Wang B S Wang J S Park H S 2015 J. Phys.: Condens. Matter 27 083001
[8] Lee C Yan H Brus L E Heinz T F Hone J Ryu S 2010 ACS Nano 4 2695
[9] Gorbachev R V Riaz I Nair R R Jalil R Britnell L Belle B D Hill E W Novoselov K S Watanabe K Taniguchi T et al. 2011 Small 7 465
[10] Molina-Sanchez A Wirtz L 2011 Phys. Rev. 84 155413
[11] Song Q J Tan Q H Zhang X Wu J B Sheng B W Wan Y Wang X Q Dai L Tan P H 2016 Phys. Rev. 93 115409
[12] Giannozzi P Baroni S Bonini N Calandra M Car R Cavazzoni C Ceresoli D Chiarotti G L Cococcioni M Dabo I Dal Corso A et al. 2009 J. Phys.: Condens. Matter 21 395502
[13] Perdew J P Zunger A 1981 Phys. Rev. 23 5048
[14] Hartwigsen C Goedecker S Hutter J 1998 Phys. Rev. 58 3641
[15] Monkhorst H J Pack J D 1976 Phys. Rev. 13 5188
[16] Hestenes M R Stiefel E 1952 J. Res. Natl. Bur. Stand. 49 409
[17] Baroni S de Gironcoli S Dal Corso A Giannozzi P 2001 Rev. Mod. Phys. 73 515
[18] Lazzeri M Mauri F 2003 Phys. Rev. Lett. 90 036401
[19] Gonze X Charlier J C Allan D C Teter M P 1994 Phys. Rev. 50 13035
[20] Gonze X Lee C 1997 Phys. Rev. 55 10355
[21] Becker D Beck H P 2004 Zeitschrift fr anorganische und allgemeine Chemie 630 41
[22] Jones R E Templeton D H 1955 Acta Cryst 8 847
[23] Wang J Dong B J Guo H H Yang T Zhu Z Hu G Saito R Zhang Z D 2017 Phys. Rev. 95 045404
[24] Nair R R Blake P Grigorenko A N et al. 2008 Science 320 1308
[25] Merthe D J Kresin V V 2016 Phys. Rev. 94 205439
[26] Tan P H Han W P Zhao W J Wu Z H et al. 2012 Nat. Mater. 11 294
[27] Stenger I Schué L Boukhicha M Berini B Plaçis B Loiseau A Barjon J 2017 2D Mater. 4 031003